%testCase for the Corpus callosum shape
headerName = 'cc.';
tailName = '.lpts';
N = 32;
n = 64;
path = 'cc-shapes/';

%store all shape in a 3D matrix===========================================
allShape = zeros(2,n,N);

for i=0:(N-1)
    filename = strcat(headerName, strcat(sprintf('%02d',i),tailName));
    allShape(:,:,i+1) = loadShapeFromFile(strcat(path,filename));
end

for i=0:(N-1)
    filename = strcat(headerName, strcat(sprintf('%02d',i),tailName));
    allShape(:,:,i+1) = loadShapeFromFile(strcat(path,filename));
end

%plot all point in all shape==============================================
figure;plotAllShape(allShape, N);

%preshape is the normalized shape without rotation aligment
%get preShape=============================================================
allpreShape = zeros(2,n,N);
for i=1:N
    allpreShape(:,:,i) = getPreShape(allShape(:,:,i));
end
figure;plotAllShape(allpreShape, N);

%OPA======================================================================
allOPAShape = zeros(2,n,N);
refShape = allpreShape(:,:,1);
for i=1:N
    allOPAShape(:,:,i) = rotateToRefShape(allpreShape(:,:,i), refShape);
end
figure;plotAllShape(allOPAShape, N);

%GPA======================================================================
allGPAShape = allpreShape;%zeros(2,n,N);
%initial estimate shape
estShape = allpreShape(:,:,1);
estShapeOld = zeros(2,n);
while(sum(sum((estShape - estShapeOld).^2))>0.001)
    for i=1:N
        allGPAShape(:,:,i) = rotateToRefShape(allGPAShape(:,:,i), estShape);
    end
    estShapeOld = estShape;
    %Procruste mean shape
    estShape = getProcrusteMeanShape(allGPAShape);
    disp('one pass finished')
end
figure;plotAllShape(allGPAShape, N);
figure; scatter(estShape(1,:),estShape(2,:),'.');

%PCA for shape estimation=================================================
D = zeros(n*2,N);
for i=1:N
    D(:,i) = [allGPAShape(1,:,i)-estShape(1,:) allGPAShape(2,:,i)-estShape(1,:)]';
end
S_trick = 1/(N-1)*(D'*D);
[V_trick, lambda] = eig(S_trick);
sqrt_lambda = sqrt(diag(lambda));
V = zeros(2*n,N);
for i=1:N
    V(:,i) = D*V_trick(:,i)/norm(D*V_trick(:,i));
end

%shape reconstruct, meanshape is the estShape
Vl=V*sqrt_lambda;

L = 0.1;
shapeL1=estShape;
shapeL1(1,:) = estShape(1,:)+ L*Vl(1:n)';
shapeL1(2,:) = estShape(2,:)+ L*Vl(n+1:2*n)';
figure; 
hold on;
scatter(estShape(1,:),estShape(2,:),'.');
scatter(shapeL1(1,:),shapeL1(2,:),'.');
hold off;
